classdef ClassifierChainsUsingLabelPriorAndPosteriorPower ...
        < MultiLabelAlgorithm ...
        & LabelSequenceable ...
        & MultiLabelAlgorithm0able ...
        & ZAlgorithmable ...
        & DecisionThresholdable ...
        & MultiLabelMetricable
    %CLASSIFIERCHAINSUSINGLABELPRIORANDPOSTERIORPOWER Summary of this class goes here
    %   Detailed explanation goes here
    
    properties ( SetAccess = protected )
        zWeight = [1/3, 1/3, 1/3] % Weight for metric in step alpha, beta and gamma
    end
    
    methods
        function [  ] = setZWeight( this, zWeight )
            if nargin == 1
                this.zWeight = [1/3, 1/3, 1/3];
            end
            if nargin == 2
                this.zWeight = zWeight;
            end
        end
    end
    
    methods
        function [ this ] = ...
                ClassifierChainsUsingLabelPriorAndPosteriorPower( name )
            if nargin == 0
                this.setName('ccul3p');
            end
            if nargin == 1
                this.setName(name);
            end
        end
    end
    
    methods
        %%
        function [  ] = build( this, X1, Y1 )
            
            %% Decide label sequence
            randomIndex = randperm(size(X1, 1));
            X = X1(randomIndex, :);
            Y = Y1(randomIndex, :);
            [nObs, nLab] = size(Y1);
            ttPercent = 0.5;
            nTt = round(nObs*ttPercent);
            
            lucMat = calcLabelUncertaintyMatrix(Y);
            lucMat = lucMat + diag(nan(size(lucMat, 1), 1));
            
            %%
            labSeq = zeros(1, nLab);
            for iLoop = 1:(nLab - 1); %% or $(nLab - 1)?
                %     keyboard %%
                metricAlpha_i = nan(size(labSeq));
                metricBeta_i = nan(size(labSeq));
                metricGamma_i = nan(size(labSeq));
                metricGamma0_i = nan(size(labSeq));
                metricGamma1_i = nan(size(labSeq));
                for jLab = find(~labSeq) ...
                        % consider the j-th label which has not been selected yet
                    % Step Alpha
                    % {X, Y_known} -predict-> y_i
                    % Generate tt(tv)Alpha_ij
                    [ttXAlpha_ij, ttYAlpha_ij, tvXAlpha_ij, tvYAlpha_ij] = ...
                        genTtTvAlpha_ij(X, Y, labSeq, jLab, nTt, lucMat);
                    algorithm_ij = this.algorithm0.algorithm0.clone();
%                     keyboard%%
%                     algorithm_ij = this.algorithm0.algorithm0.clone();
                    algorithm_ij.build(ttXAlpha_ij, ttYAlpha_ij);
                    result_ij = algorithm_ij.apply(tvXAlpha_ij, tvYAlpha_ij);
%                     keyboard%%
                    metricAlpha_i(jLab) = this.metric.apply(tvYAlpha_ij, result_ij.Y_hat, result_ij.Y_out);
%                     modelAlpha_ij = feval(hSubTrainFcn, ...
%                         ttXAlpha_ij, ttYAlpha_ij, subTrainOption);
%                     resultAlpha_ij = feval(hSubApplyFcn, ...
%                         modelAlpha_ij, tvXAlpha_ij, tvYAlpha_ij, subApplyOption);
%                     tvYhatAlpha_ij = resultAlpha_ij.Yhat;
%                     metricAlpha_i(jLab) = calcHammingLoss(tvYAlpha_ij, tvYhatAlpha_ij);
                    
                    % step Beta
                    % {X, Y_known, y_i} -predict-> Y_unknown
                    % Genearate tt(tv)Beta_ij.
                    [ttXBeta_ij, ttYBeta_ij, tvXBeta_ij, tvYBeta_ij] = ...
                        genTtTvBeta_ij(X, Y, labSeq, jLab, nTt, lucMat);
                    % do validation Beta
                    algorithm_ij = this.algorithm0.clone();
                    algorithm_ij.build(ttXBeta_ij, ttYBeta_ij);
                    result_ij = algorithm_ij.apply(tvXBeta_ij, tvYBeta_ij);
                    metricBeta_i(jLab) = this.metric.apply(tvYBeta_ij, result_ij.Y_hat, result_ij.Y_out);
%                     modelBeta_ij = feval(hSubTrainFcn, ...
%                         ttXBeta_ij, ttYBeta_ij, subTrainOption);
%                     resultBeta_ij = feval(hSubApplyFcn, ...
%                         modelBeta_ij, tvXBeta_ij, tvYBeta_ij, subApplyOption);
%                     tvYhatBeta_ij = resultBeta_ij.Yhat;
%                     metricBeta_i(jLab) = calcHammingLoss(tvYBeta_ij, tvYhatBeta_ij);
                    
                    % step Gamma
                    % Suppose that the j-th label is unknown here.
                    % Check whether is the j-th label dependent is on the remaining
                    % other labels?
                    % Gamma0
                    metricGamma0_i(jLab) = metricAlpha_i(jLab);
                    % Gamma1
                    [ttXGamma1_ij, ttYGamma1_ij, tvXGamma1_ij, tvYGamma1_ij] = ...
                        genTtTvGamma1_ij(X, Y, labSeq, jLab, nTt, lucMat);
                    % do validation Gamma1
                    algorithm_ij = this.algorithm0.algorithm0.clone();
                    algorithm_ij.build(ttXGamma1_ij, ttYGamma1_ij);
                    result_ij = algorithm_ij.apply(tvXGamma1_ij, tvYGamma1_ij);
                    metricGamma1_i(jLab) = this.metric.apply(tvYGamma1_ij, result_ij.Y_hat, result_ij.Y_out);
%                     modelGamma1_ij = feval(hSubTrainFcn, ...
%                         ttXGamma1_ij, ttYGamma1_ij, subTrainOption);
%                     resultGamma1_ij = feval(hSubApplyFcn, ...
%                         modelGamma1_ij, tvXGamma1_ij, tvYGamma1_ij, subApplyOption);
%                     tvYhatGamma1_ij = resultGamma1_ij.Yhat;
%                     metricGamma1_i(jLab) = calcHammingLoss(tvYGamma1_ij, tvYhatGamma1_ij);
                    
                end
                % Gamma0 - Gamma1 = unknown - known
                %               = the help of knowing the remaining labels
                % The smaller this difference is, the less j-th label is dependent
                % on the remaining labels, i.e., the better
                metricGamma_i = abs(metricGamma1_i - metricGamma0_i);
                
                [~, sortedAlphaId_i] = sort(metricAlpha_i, 'descend');
                [~, sortedBetaId_i] = sort(metricBeta_i, 'descend');
                [~, sortedGammaId_i] = sort(metricGamma_i, 'ascend');
                sortedBetaId_i(isnan(metricBeta_i)) = -inf;
                sortedGammaId_i(isnan(metricBeta_i)) = -inf;
%                 keyboard%%
                [~, bestId_i] = min( ...
                    (this.zWeight(1)*sortedAlphaId_i + ...
                    this.zWeight(2)*sortedBetaId_i + ...
                    this.zWeight(3)*sortedGammaId_i)/sum(this.zWeight));
                
                %     metricGamma_i = metricGamma1_i - metricGamma0_i;
                %     [~, sortedAlphaId_i] = sort(metricBeta_i - metricAlpha_i, 'ascend');
                %     [~, sortedBetaId_i] = sort(metricBeta_i, 'ascend');
                %     [~, sortedGammaId_i] = sort(metricGamma_i, 'ascend');
                %     sortedBetaId_i(isnan(metricBeta_i)) = inf;
                %     sortedGammaId_i(isnan(metricBeta_i)) = inf;
                %     [~, bestId_i] = min( ...
                %         (weight(1)*sortedAlphaId_i + ...
                %         weight(3)*sortedGammaId_i)/(weight(1) + weight(3)));
                
                labSeq(bestId_i) = iLoop;
            end
            labSeq(labSeq == 0) = nLab;
            
            this.labelSequence = labSeq;
            
            %% Build sub-algorithms
            
            X1_augment_i = X1;
            for iLabel = this.labelSequence
                Y1_i = Y1(:, iLabel);
                
                this.zAlgorithm{iLabel} = this.algorithm0.clone();
                this.zAlgorithm{iLabel}.build(X1_augment_i, Y1_i);
                
                X1_augment_i = [X1_augment_i, Y1_i];
            end
        end
        
        %%
        function [ result ] = apply( this, X2, Y2 )
            
%             nLabel = size(Y2, 2);
            
            result.Y_hat = zeros(size(Y2));
            result.Y_out = result.Y_hat;
            
            X2_augment_i = X2;
            for iLabel = this.labelSequence
                Y2_i = Y2(:, iLabel);
                
                result_i = this.zAlgorithm{iLabel}.apply(X2_augment_i, Y2_i);
                
                result.Y_hat(:, iLabel) = result_i.Y_hat;
                result.Y_out(:, iLabel) = result_i.Y_out;
                
                X2_augment_i = [X2_augment_i, result_i.Y_hat];
            end
        end
    end
    
end

%%
function [ liecMat ] = calcLIECMatrix( Y )
% Calculate Label Inter-Effect Coefficient Matrix

[nObs, nLab] = size(Y);

liecMat = zeros(nLab, nLab);
for iLab = 1:nLab
    for jLab = 1:nLab
        y_i = Y(:, iLab);
        y_j = Y(:, jLab);
        id0 = y_i == 1;
        %         id0 = ~(y_i == -1 & y_j == -1);
        id1 = (y_i == 1 & y_j == 1);
        liec_ij = nnz(id1)/nnz(id0);
        if isnan(liec_ij)
            liec_ij = -1;
        end
        liecMat(iLab, jLab) = liec_ij;
    end
end

end

%%
function [ Z ] = calcLabelUncertaintyMatrix( Y )

Z = calcLIECMatrix(Y);

Z1 = Z.*log2(Z);
Z2 = (1 - Z).*log2(1 - Z);

Z = -(Z1 + Z2);
Z(isnan(Z)) = 0;

end

%%
function [ttXAlpha_ij, ttYAlpha_ij, tvXAlpha_ij, tvYAlpha_ij] = ...
    genTtTvAlpha_ij(X, Y, labSeq, jLab, nTt, lucMat)
% lucMat - label uncertainty matrix

if jLab; end
nLab = length(labSeq);
knownYAlphaId_ij = false(1, nLab);
knownYAlphaId_ij(logical(labSeq)) = true;
unknownYAlphaId_ij = ~knownYAlphaId_ij;

% % Select only a subset of Y_known which may help most to predict y_i
% luc_j = lucMat(:, jLab);
% knownYAlphaId_ij = knownYAlphaId_ij & (luc_j <= nanmean(lucMat(:)))';

% training training Alpha
ttXAlpha_ij = horzcat( ...
    X(1:nTt, :), ...
    Y(1:nTt, knownYAlphaId_ij));
% ttYAlpha_ij = Y(1:nTt, unknownYAlphaId_ij);
ttYAlpha_ij = Y(1:nTt, jLab);

% training validation Alpha1
tvXAlpha_ij = horzcat( ...
    X((nTt + 1):end, :), ...
    Y((nTt + 1):end, knownYAlphaId_ij));
% tvYAlpha_ij = Y((nTt + 1):end, unknownYAlphaId_ij);
tvYAlpha_ij = Y((nTt + 1):end, jLab);

end

%%
function [ ttXBeta_ij, ttYBeta_ij, tvXBeta_ij, tvYBeta_ij ] = ...
    genTtTvBeta_ij( X, Y, labSeq, jLab, nTt, lucMat )

% suppose that the j-th label is known here
knownYBetaId_ij = logical(labSeq);
knownYBetaId_ij(jLab) = true; % assume that the j-th label is known here
% knownYBetaId_ij = logical(knownYBetaId_ij);
unknownYBetaId_ij = ~knownYBetaId_ij;

luc_j = lucMat(:, unknownYBetaId_ij);
knownYBetaId_ij = knownYBetaId_ij & (any(luc_j < nanmean(lucMat(:)), 2))';

% training training Beta
ttXBeta_ij = horzcat( ...
    X(1:nTt, :), ...
    Y(1:nTt, knownYBetaId_ij));
ttYBeta_ij = Y(1:nTt, unknownYBetaId_ij);
% training validation Beta
tvXBeta_ij = horzcat( ...
    X((nTt + 1):end, :), ...
    Y((nTt + 1):end, knownYBetaId_ij));
tvYBeta_ij = Y((nTt + 1):end, unknownYBetaId_ij);

end

%%
% function [ ttXGamma0_ij, ttYGamma0_ij, tvXGamma0_ij, tvYGamma0_ij ] = ...
%     genTtTvGamma0_ij( X, Y, labSeq, jLab, nTt )
%
% if jLab; end
% nLab = length(labSeq);
% knownYGammaId0_ij = false(1, nLab);
% knownYGammaId0_ij(logical(labSeq)) = true;
% unknownYGammaId0_ij = ~knownYGammaId0_ij;
% % training training Gamma1
% ttXGamma0_ij = horzcat( ...
%     X(1:nTt, :), ...
%     Y(1:nTt, knownYGammaId0_ij));
% ttYGamma0_ij = Y(1:nTt, unknownYGammaId0_ij);
% % training validation Gamma1
% tvXGamma0_ij = horzcat( ...
%     X((nTt + 1):end, :), ...
%     Y((nTt + 1):end, knownYGammaId0_ij));
% tvYGamma0_ij = Y((nTt + 1):end, unknownYGammaId0_ij);
%
% end

%%
function [ ttXGamma1_ij, ttYGamma1_ij, tvXGamma1_ij, tvYGamma1_ij ] = ...
    genTtTvGamma1_ij( X, Y, labSeq, jLab, nTt, lucMat )

nLab = length(labSeq);
knownYGammaId1_ij = true(1, nLab);
knownYGammaId1_ij(jLab) = false;
unknownYGammaId1_ij = ~knownYGammaId1_ij;

luc_j = lucMat(:, unknownYGammaId1_ij);
knownYGammaId1_ij = knownYGammaId1_ij & (luc_j < nanmean(lucMat(:)))';

% training training Gamma1
ttXGamma1_ij = horzcat( ...
    X(1:nTt, :), ...
    Y(1:nTt, knownYGammaId1_ij));
ttYGamma1_ij = Y(1:nTt, unknownYGammaId1_ij);
% training validation Gamma1
tvXGamma1_ij = horzcat( ...
    X((nTt + 1):end, :), ...
    Y((nTt + 1):end, knownYGammaId1_ij));
tvYGamma1_ij = Y((nTt + 1):end, unknownYGammaId1_ij);

end
